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I. INTRODUCTION 



The exploration of the thermodynamic behavior of materials under extreme conditions 
usually follows two paths corresponding to two experimental devices: isothermal compres- 
sions and shock compressions. Isothermal compressions are performed with diamond anvil 
cell (DAC) techniques, and are used to compress materials up to very high pressures, al- 
though with limited temperatures. On the other hand, shock compression experiments 
investigate the high-pressure/high temperature regions through the propagation of dynamic 
shock waves in the system. Nevertheless, the thermodynamic domain available using shock 
experiments remains limited to the so-called Hugoniot curve, which is, by definition, the 
collection of thermodynamic states which can be reached from a system at fixed initial con- 
ditions, with shocks of increasing strengths. Such a constraint arises from the fact that the 
system is assumed to satisfy Euler's fluid equations {i.e inviscid Navier-Stokes equations), 
and physically meaningful shocks therefore have to fulfill the Rankine-Hugoniot conditions, 
which relate the thermodynamic parameters of the fluid at rest, and the thermodynamic 
parameters of the shocked material. Another constraint is that shock waves are adiabatic, 
therefore leading to very large temperature increases in the material, which limits its com- 
pressibility. The equations of state (EOS) used to predict the material's behavior at the 
extreme conditions encountered are often simple extrapolations of EOS fitted on available 
data, i.e. shock data and DAC data. It then appeared interesting to enlarge the experimen- 
tal domain of investigation of materials behavior using dynamic compression set-ups, and 
particularly isentropic compressions. Several experimental set-ups allow to load a pressure 
ramp in a material. The first one is the high pulsed power (of which the sandia Z machine 
and the High Explosive Pulsed Power-i^ at LANL are good examples). The second one con- 
sists in using an impactor with a varying density along one direction, as proposed initially at 
the AIP - SWCM conference.- A successful technique is to stack slices of different materials, 
leading to the so-called PILLOW impactors at Sandia,- MIVAR impactors in France,- and 
more recently the FGM (Functionaly Graded Materials) impactors at LLNL,- allowing a 
real design of a thermodynamic path as a succession of shock and release waves. The last 
one concerns experiments of Barnes' type where the compression is the consequence of the 
isentropic release of another material, as for example detonation products.-i^ 

Experiments involving isentropic compression are of great interest to reach high compres- 
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sion states, or in geophysic applications to reach states representative of the earth's core. 
Experiments involving a precise evaluation of release waves in materials need also general 
numerical methods to compute the states reached by isentropic pocesses. 

Up to now, simulations involving shock processes are rather well developed due to the 
simplicity of the Hugoniot equat ions. Those studies are performed within the framework of 

for reference textbooks on molecular simulation, and Ref. 



statistical physics, see Refs.|9 



10 



11 



for reference works on nonequilibrium simulation of shock waves. Any state lying on the 
Hugoniot can be reached from the reference state by searching for a given compression the 
temperature for which the pressure and the total energy of the system satisfies the Hugoniot 
relation. The search can be implemented in very efficient manners.— >i^>i^ Those methods 
are now adapted to classical molecular dynamics and Monte Carlo, as well as quantum 
molecular dynamics.— 

Such an easy method does not exist for isentropic processes. In this paper, we present or 
recall several equilibrium methods which allow to follow isentropic paths, both for classical 
or quantum molecular dynamics simulations. We contrast these methods in terms of their 
precisions, rigor and computational requirements. We compare the results obtained from 
equilibrium simulations with release waves observed in nonequilibrium molecular dynamics. 
The comparison between equilibrium and nonequilibrium methods therefore measures how 
isentropic the expansion of the system is. It is expected that release waves of a perfect 
non-viscous fluid are isentropic. For simple monoatomic fluids such as argon, it is often 
assumed that the release is isentropic, and viscosity effects are neglected. Our results show 
that even in this simple case, the release is not strictly isentropic and some corrections have 
to be taken into account. As a by-product of our study, we also explore more precisely the 
relationship between the Hugoniot and the isentrope curves, from a numercial viewpoint, but 
also giving a statistical physics proof of the coincidence of the curves for small compressions 
(see Appendix B). 



Organization of the paper 

The paper is organized as follows. In Section [Tll we present the nonequilibrium method 
used to simulate rarefaction waves, while some equilibrium methods for constant entropy 
sampling are recalled in Section IIIII - details of the practical implementaion of the thermo- 
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dynamic integration and precisions on the isentropic integration are given in Appendix A 
and C respectively). In Section HVl a comparison of numerical results obtained in the case 
of release waves in argon is performed, and we discuss whether release waves are isentropic. 



II. NONEQUILIBRIUM SIMULATIONS 

A. Microscopic description of physical systems 

We first recall how a system is described at the microscopic level using statistical 
mechanics.— Statistical physics is a theory which allows to compute thermodynamic (macro- 
scopic) properties of a system knowing the microscopic interactions between its constitutive 
elements. 

Consider a microscopic system composed of particles, confined in a simulation box 
T) = [0, Lx] X [0, Ly] X [0, Lz]. The volume of the domain is denoted hj V = \V\ = L^LyL^. 
The system is characterized by the positions q = (gi, . . . , q^) and momenta p = (pi, . . . ,pn) 
of the particles, which have masses mj, and interact through a potential energy function U. 
The phase-space Q is the collection of all possible microscopic configurations {q,p) of the 
system. 

The central quantity describing the system is the Hamiltonian 



which gives the energy of a given microscopic configuration {q,p). Average thermodynamic 
properties of the system can be computed as averages of functions of the microscopic vari- 
ables 0{q,p) (the so-called observables) with respect to the canonical measure at a temper- 
ature T, for a given simulation box: 



The canonical measure associated with the Hamiltonian ([T]) weights microscopic states ac- 
cording to their energies using a Boltzmann weight:— 
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where /cb is the Boltzmann constant, and the partition function Zy^T is a normahzation 
factor so that ([3]) is indeed a probabihty measure: 



We have indicated exphcitly the dependence of the canonical measure ([3]) and the partition 
function (j4]) on the temperature T and the volume V since these will be the parameters 
allowed to vary in the sequel. 

B. Nonequilibrium simulation of release waves 

Similarly to what has been proposed for the simulation of shock waves,— isentropic com- 
pressions or releases can be simulated directly using Non-Equilibrium Molecular Dynamics 
(NEMD). A straightforward numerical set-up to this end is simply to throw a low speed pis- 
ton towards the sample (creating a weak shock), and then accelerating the piston in time. 
Except this external forcing, the system evolves according to the standard hamiltonian dy- 
namics 



which is integrated in time with the Verlet scheme.— A linear compression ramp would be 
obtained in the case where the acceleration is constant in time. 

To obtain release waves, a shock wave can be loaded in a sample; when this shock wave is 
reflected when interacting with a free surface, it transforms into a release wave, supposedly 
isentropic. In this study, we start the release from an equilibrated state obtained from 
a preliminary canonical simulation, using three dimensional periodic boundary conditions. 
When the system is equilibrated, the periodic boundary conditions are removed in the x 
direction. Two release waves are then created at the two free surfaces, and they propagate 
in opposite directions towards the center of the box. This process is illustrated in Figure [H 

From the simulation data presented in Figure [T], profiles of thermodynamic quantities 
(average densities, (kinetic) temperatures and pressures) can be extracted and averaged 
over thin slices. Moreover, the two release waves being symmetric, their related profiles can 
be averaged. A superposition of the profiles, taken at different times but projected back in 




(4) 




(5) 
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the same thermodynamic diagram, is then obtained and averaged over, leading to a single 
profile. 

The accuracy of the computation increases with the system size: an increase in the size 
of the transverse directions decreases the uncertainties on the slice averages (thanks to a 
thermodynamic limit), and an increase in the longitudinal direction allows to accumulate 
more profiles in time, therefore reducing statistical errors. 

III. EQUILIBRIUM METHODS FOR CONSTANT ENTROPY SAMPLING 

We present in this section three methods to compute the collection of all states (in terms 
of their temperature, volume and pressure) which have the same entropy as some reference 
state. These methods therefore allow to draw a curve in the (V, T) diagram (or in the 
(P, T) or {P,V) diagrams), called the isentrope, and will be used as benchmark methods 
in Section [IV] to check whether release waves computed by NEMD simulations are indeed 
isentropic or not. We emphasize that, altough presented for the computation of isentropic 
releases, all the methods described in this section may also be used to determine isentropic 
compressions. We also recall a fourth method, used to obtain the entropy of a system once 
the entropy of some reference state (such as the perfect gas) is fixed. The computational 
cost of the latter method as well as its low accuracy for dense states prevented us from 
applying it to enough points to obtain an entire isentrope curve, and we therefore limited 
its use to a consistency check on the results obtained with the other methods. 

A. Thermodynamic integration 

The entropy of the system varies when the simulation conditions are changed. Here, we 
consider that the states visited by the release wave are a succession of local equilibrium 
states, which can be described within the canonical ensemble as given by statistical physics. 
Therefore, the state of the system is defined by two parameters, its volume (equivalently, 
the density) and its temperature. 
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1. Variables indexing the variations. 



Consider a general transformation in which both the volume accessible to the system 
(equivalently, the density) and the temperature are varied. We restrict ourselves to variations 
of the domain in one spatial direction only, to model the anisotropic behavior of release waves. 
Assuming that the state of the system at rest can be described by some cubic simulation 
box with periodic boundary conditions, the volume under compression may be indexed by 
a variable Ai, so that the associated simulation domain 'D{\i) = [0, (1 + Xi)Lx] x [0, L]^ has 
a volume 

V{X,) = {1 + X,)L,L\ 

Notice that we consider ^ L since we may start from a uniaxially compressed state. The 
temperature variations are indexed by a parameter A2: 

T(A2) = (1 + A2(5T)T, 

for some reference temperature T and a given relative temperature variation ST, the tem- 
perature variation being therefore AT = T6T. The reference inverse temperature is still 
(3~^ = k-QT. The particular case where only the temperature is changed (while the volume 
is kept constant) corresponds to Ai constant, while isothermal transformations are charac- 
terized by A2 remaining constant. Expansions correspond to Ai > 0. 



2. Parametrization of the isentrope curve. 

The isentrope is the locus of the points in the (Ai, A2) space such that the entropy nor- 
malized by the Boltzmann factor 

^ = ^ (6) 

ks ksT 

is constant, JF denoting the free energy of the system, and U its energy. This thermodynamic 
relation can be converted into an equivalent formula in the framework of statistical physics, 
which is much more convenient from a computational viewpoint: 

l(=V({T, V) = {H)v,T, 

where the canonical average is defined in Eq. ([2]), and 



Jn 



■(3Hiq,p) 



dqdp. 



We start from some reference state described by the parameters (Ai,A2) = (0,0). The 
statistical physics reformulation of the requirement that ([6]) be constant is then 



\ Jq{o) / 
In this expression, the phase-space fi(Ai) is the collection of all possible microscopic config- 
urations of the system associated with a domain T>{Xi) of volume V{Xi). 

3. Numerical implementation 

To determine the isentrope curve, we compute the entropy variation along a given path 
in (Ai, A2) space going through the reference initial state, and search for the point such that 
the entropy difference with this state is 0. A simple choice is illustrated in Fig. [2l It consists 
in performing 

(i) an isothermal rarefaction, going from the initial compressed state (0, 0) to an interme- 
diate state (Ai,0) with Ai > 0; 

(ii) in a second step, an isochore cooling, going from the intermediate state (Ai, 0) to some 
final state (Ai,A2), resorting to a maximal temperature difference A2AT < large 



The idea is that, in general, the first part of the transformation increases the entropy of the 
system (since more space becomes available for the particles), while the entropy decreases 
in the second part (since the temperature decreases). Of course, more general paths, with 
joint variations of Ai and A2, could be considered. 

The energies {H)y(^Xi),t{X2) ^t^^ computed using standard sampling strategies, while the 
remainder in the expression of S{Xi, A2) —5(0, 0), a ratio of partition functions, is estimated 
using standard techniques for free-energy calculations. This is detailed in Appendix A. 

We emphasize that this procedure is time consuming since it requires many equilibrium 
samplings to obtain one point on the curve. It is however exact (up to statistical errors and 
discretization errors in the integrals defining A), and can be straightforwardly parallelized 
since the equilibrium samplings required are independent. 



5(Ai,A2)- 5(0,0) 




{H)v(0),T{0) 



enough. 
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B. Successive hugoniostat simulations 



The variations of macroscopic quantities across a shock interface are governed by the 
Rankine-Hugoniot relations, which relate the jumps of the quantities under investigation 
(pressure, density, velocities) to the velocity of the shock front. The third Rankine-Hugoniot 
conservation law for the Euler equation governing the hydrodynamic evolution of the fluid 
reads (macroscopic quantities are denoted by curly letters) 

n = u-u^-]^{v + Vo){Vo-v) = 0. (7) 

In this expression, U is the internal energy of the fluid, V its pressure, and V its volume. 
The subscript refers to the initial state (the pole), the other quantities are evaluated at 
a state obtained from some shock compression, after equilibration. The Hugoniot curve 
corresponds to all the possible states satisfying ([7]). In practice, the collection of these states 
may be computed by nonequilibrium simulations with shocks of different strengths, inducing 
various compressions. 

Alternatively, small equilibrium simulations may be used, relying on the statistical physics 
reformulation of the Hugoniot relation: 

7Y(Ai, A2) - 7Y(0, 0) = (-f/')\/(Ai),T{A2) - (-f^)v'(0),T(0) 

+ Y^(0) ((P..)y(AO,T(A.) + (P..)y(0),T(0)) = 0. (8) 

The XX component of the pressure tensor is, for a simulation domain of volume ^(Ai), 

1 ^ 
^(^1) 

For a given variation of the volume for instance (indexed by Ai), the variation \2T6T of the 
temperature is sought for, using for instance the techniques described in Refs. Il2|jl3 . 

The Hugoniot curve does not have a priori any relationship with the isentrope curve. 
However, it can be shown that the entropy variation along the Hugoniot curve is negligible 
up to terms of order three in the volume variable; the Hugoniot and the isentropic curves 
are osculatory. We present in Appendix B two proofs, the standard proof based on thermo- 
dynamic relations, and a new proof fully relying on a statistical physics reformulation. 

The good agreement between the Hugoniot and the isentrope for small compressions 
and/or expansions can be used to compute the isentropic curve as a succession of weak shocks 



or weak releases, this approximation getting more accurate as the shock compressions are 
weakened. The only parameter left in this method is the relative volume change W jV = 
A""*"^ — A" during the instantaneous compressions or releases. We used the Hugoniostat 
methodi^ii^ to compute a sequence of states (A", A2) such that 7Y(A"^\ -^2^^) = ^(-^i, -^2)! 
the corresponding thermodynamic properties at these states being obtained as a by-product 
of the simulation. 



C. Isentropic integration 



Another way to perform thermodynamic integration along an isentropic path has been 

dV 

proposed by Desjarlais.— The method relies on the equilibrium evaluation of (see 
Eq. f|TT]) below). It could be applied to a system where the pressure is not isotropic upon 
replacing the pressure observable by the xx component of the pressure tensor. 
The total differential of the entropy can be written as: 



"5 = 1^ 

or 



dT+ — 
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For constant entropy processes, 
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so that, along the isentrope, 

dT _ 
f ~ 

This equation can be integrated as 



dV 




df 
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dU 
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V 



dV. 



T2 

— = exp 



V2 



Vi 



dV 

du 



dV 



(10) 



(11) 



giving the temperature T2 at which the system at volume V2 has the same entropy as the 
system in the reference state {Ti^Vi). This formula is evaluated in practice by discretizing 
the integral appearing in the exponential, and approximating the integrand using standard 
canonical sampling procedues. We refer to Appendix C for more precisions. 
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D. Evaluation of the entropy based on the chemical potential 



This technique, which can be used only for systems in a fluid phase, follows the classical 
methodology of computing the free energy JF of a system starting from the thermodynamic 
relation:— 

= U -TS = N^i-VV, (12) 
where /i is the chemical potential, defined in the canonical ensemble as 

In the case of a canonical simulation, all thermodynamic quantities are functions of the 
volume and the temperature, so that 

^ ^ U{T,V) + N^^{T,V)-V{T,V)V ^ ^^^^ 

This expression allows to compute the absolute entropy of the system provided the chemical 
potential is known,— the average pressure and energy being computed using standard sam- 
pling techniques. The chemical potential is estimated using the Widom insertion method. 



IV. NUMERICAL RESULTS FOR RELEASE WAVES 

We compare in this section the results for the different techniques presented in Sections IXTl 
and mil for a release in a Lennard- Jones system (argon). The aim is to assess whether the 
release is indeed isentropic, and also to demonstrate that approximate equilibrium compu- 
tations for small systems (successive Hugoniostat, isentropic integration) can approximate 
the isentrope curve obtained from the more rigourous and costly thermodynamic integration 
technique. 



A. Numerical parameters 



1. Initial state. 



We consider argon in an initial shocked state, located on the Hugoniot curve for a 
compression such that = cL with c = 0.65, and corresponding to T = 1758 K and 
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P = 1.7 X 10^'^ Pa. At these thermodynamic conditions, the system is in a hquid state. The 
interactions within noble gas atoms are well- described by a Lennard- Jones potential: 

l<i<j<N ^ ^ 

In the case of Argon, e/ks = 120 K and a = 3.405 A. The cut-off radius for the Lennard- 
Jones interaction is here r^ut = 2.5a. 

2. Nonequilihrium simulations. 

In order to reach this initial state before performing the NEMD release, a preliminary 
hugoniostat simulation is run for a system of 50 x 50 x 500 unit cells, using periodic boundary 
conditions. Then, the boundary conditions in the longitudinal direction are removed, and 
the system evolves according to the Hamiltonian dynamics. Profiles of thermodynamic 
quantities are computed every 0.25 ps for the post-processing procedure described at the 
end of Section III B[ 

3. Equilibrium simulations 

Equilibrium computations have been performed with a system composed of = 4000 
atoms, starting in a FCC crystal geometry before melting, using periodic boundary condi- 
tions in all directions. 

a. Thermodynamic integration. As shown in Section IIII At the search of states having 
the same entropy can be performed using thermodynamic integration, which amounts to per- 
forming many equilibrium simulations. The canonical sampling for a given set of parameters 
(Ai, A2) is done with a Langevin dynamics for Agtcps = 2^^ time steps, with At = 2 x 10^^^ s, 
and a friction coefficient 7 = 10^'^ s~^. 

First, the entropy variation along the isothermal release is computed, with canonical 
samplings along the path (0, 0) (0, Ai) with Ai = 0.54 (using M + 1 = 15 states). Then, 
for each compression of interest, the isochore cooling is performed using temperature steps 
AT = —25 K for expansions Ai < 0.25, and AT = —50 K for states Ai > 0.25 (these paths 
can be restated in terms of A2 G [0, 1] upon considering a temperature modification AT 
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depending on the compression). The numerical integration for computing the value of A is 
finally performed using the trapezoidal rule. 

Error estimates on the canonical samplings are obtained with block averaging.— In all 
the cases considered, the statistical error (as measured using the 95% confidence interval 
associated with the variance computed from block-averaging) is inferior to 1%. Therefore, 
the entropy difference is computed within 1% errors. For fixed Ai, the state A2 such that 
<S{Xi, A2) is constant is then known with an error depending on the local value of the partial 
derivative of S with respect to A2. This error can immediately be reformulated as an error 
on the estimated temperature. The error on the computed pressure is the error arising from 
the error on the state A2, plus the sampling error. It is found to be at most 2%. 

b. Successive Hugoniostat. Successive Hugoniostat simulations have been performed 
with a Langevin version of the Hugoniostat method (see Eq. (11) in Ref. 



22 



with the 



parameters ^ = 10^^ s~^ and u = 10^^ s~^). Trajectories of A'steps = 50,000 timesteps at 
each compression are considered, with a timestep At = 5 x 10~^^ s. The relative volume 
change SV/Vq from one point on the curve to another is set to 0.01. 

c. Isentropic integration. See Appendix C. 

d. Entropy evaluation. The test particle insertion method used to evaluate the chem- 
ical potential requires many more iterations than the other equilibrium techniques. In the 
same framework as for isentropic integration (see Appendix C), A'^stcps = 5 x 10^ iterations 
were needed to obtain a satisfactory convergence. The statistical error on the calculated en- 
tropy (as measured using the 95% confidence interval associated with the variance computed 
from block-averaging) is estimated to be inferior to 1.2 %. 



B. Discussion of the numerical results 



Release waves are presented in Figures [3H5] in three different diagrams, {P, p),{P,T) and 
{T,p). It can clearly be seen that the results coming from the three equilibrium techniques 
of isentropic simulations are very close. This shows that provided the relative volume change 
parameter is carefully chosen in either the successive hugoniostat method or the isentropic 
integration, the propagating error remains at a low level; these methods can then be as 
accurate as the more rigorous and costly method of thermodynamic integration. Moreover, 
evaluating the chemical potential, we have computed absolute value of the entropy at three 
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densities, p = 2780 kg.m ^, p = 2190 kg.m ^, and p = 1806 kg.m ^. The corresponding 
values, 83.2±0.95 J.mor\ 83.1 ±0.42 J.mol-\ 83.2±0.18 J.mor\ confirm that the entropy 
is indeed constant (within the error bars) on the calculated curve, validating once again the 
different methods. 

The comparison with the results of the expansion of the liquid using nonequilibrium MD is 
also fruitful. The overall agreement is fair enough, which means that release waves are indeed 
almost isentropic. However, it can be noticed that the temperature is not predicted correctly. 
While the different curves look very similar in a (P, p) diagrams, some discrepancies appear 
in the (T, P) diagram, which are even more obvious in the (T, p) diagram. Indeed, for the 
latter diagram, the observed temperatures around the final density are greater than the 
error bars. The thermodynamic path followed by the system during its release exhibits 
systematically a higher temperature than the one of an isentropic process. This means that 
the release of a monoatomic liquid is not strictly isentropic, as is sometimes expected or 
assumed. 

Recall however that a release is expected to be isentropic only for non-turbulent flows 
of non-viscous fluids. In the case considered here, the fluid has a finite, non zero viscosity, 
and therefore dissipates energy under the form of heat. As a consequence, the temperature 
should be higher than for an isentropic release. A tentative of evaluation of this effect is 
presented below. On the Hugoniot curve,— viscous effects can be introduced in the Hugoniot 
relation^ by means of the "viscous pressure" tt as 

^7r(Vo -V)=U-Uo-^{V + Vo)iVo - V), (15) 

where tt is defined as 

du 

V being the fluid viscosity and — the velocity gradient. Taking the viscosity of the argon 

ax 

fluid at T = 700 K and P = 1 GPa (the most extremes conditions of available thermody- 
namic tables), and considering an average velocity gradient (taken during the fluid release), 
we find a temperature elevation of a few Kelvins. Considering that the pressure is much 
higher in our simulation, and therefore that the viscosity should be also greater, the actual 
temperature increase due to the finite viscosity should rather be of the order of a few tens 
of Kelvins, which is consistent with what can be observed in our numerical results. 
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Finally, a purely nonequilibrium effect has been observed during the NEMD simulations, 
that also leads to a temperature increase in the system. Indeed, gradients of thermodynamic 
and kinematic quantities are large at the first stages of the release, when the hot and dense 
material is in contact with void. The thermodynamic path followed by the system at those 
early stages of the simulation does not correspond to the thermodynamic path followed when 
the release has reached its self-similarity regime. Some equilibration time is needed for some 
steady-state regime to be reached. We evaluated this time to be around 5 picoseconds. 

V. CONCLUSION 

We have presented or recalled several equilibrium methods to compute isentropic pro- 
cesses in the high pressure regime, either for compressions or releases. These methods, 
although very different in nature, lead to similar results when applied to the release of a 
monoatomic liquid. 

We have then compared release waves computed with these equilibrium methods with 
the nonequilibrium simulation of the release process. The results show that the release is 
almost, but not strictly isentropic, the system's temperature being systematically greater 
than the one of the isentropic process. This is the consequence of two effects. First, the fluid 
actually has a finite viscosity and therefore dissipates heat, leading to a temperature increase. 
To our knowledge, this is the first time that this effect has been quantified rigorously using 
molecular dynamics simulations. Moreover, the thermodynamic path followed by the system 
during its release takes some time to reach a converged profile. We anticipate that these 
effects will be enhanced in the case of a more complex fluid, for example in the case of a the 
release of detonation product. Therefore, the assumption that release waves are isentropic 
should be carefully verified in each case. 
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APPENDIX A: PRACTICAL IMPLEMENTATION OF THE THERMODY- 
NAMIC INTEGRATION 



1. Reformulation of the problem in a fixed geometry 

From a computational viewpoint, it is more convenient to work with a fixed simulation 
domain. For instance, the unperturbed domain V{0) may be used to fix the geometry 
of the system. The volume variations are then rephrased as variations in the interaction 
scale between the particles in the direction of compression or release. In the same vein, 
the temperature may be kept constant, upon rescaling the interactions strengh by a factor 
depending on the temperature variation. Introducing the rescaled potential energy for a 
configuration q = (x, y, z): 



canonical averages for a volume V{Xi) at a temperature T(A2) can be reformulated as canon- 
ical averages in terms of the rescaled Hamiltonian H\^ x^ at the reference state at volume 
V^(0) and temperature T(0). More precisely. 



1 



U{{l + \i)x,y,z). 



1 + \2^T 



and the associated Hamiltonian 




{H)v{\^).T{\2) 



—kBT{X2) + {l + X2ST){{Ux,M)) 



where 




It is then easily seen that 



5(Ai,A2) -5(0,0) = — ln(l + A25T) + Arin(l + Ai) 

+/5[((f/A„A.))A„A.-((t/0,0))0,0] 



(Al) 



+^(Ai,A2), 
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with 



In 



v{oy 



V 



-/3C/o,o(9) 



dq 



v{oy 



In the above expression of the entropy difference, the first hne is the ideal gas contribution 
to the entropy difference. As a consistency check, we can verify that the entropy increases 
when the volume or the temperature is increased, as expected. The terms on the second 
and third lines in Eq. ( lAll) are the "excess" contributions associated with the potential 
interaction energy. 

2. Numerical evaluation of the difTerent terms 



To estimate S, two quantities are required: 

(i) averages ((■))ai,A2 with respect to the Hamiltonian H^^ x^ are computed using stan- 
dard sampling techniques such as a Langevin dynamics at an inverse temperature /3, 
implemented using the so-called BBK algorithm.— Of course, many other sampling 
techniques could be used to estimate this canonical average, in particular Nose-Hoover 



dynamics2ii2a or Metropolis-Hastings schemes^iiSa (see Ref. 



29 



for a mathematical re- 



view on sampling methods in the context of molecular simulation); 

(ii) the term A{\i,\2) requires more care in its estimation. Since this term is a ratio 
of partition functions, standard techniques used for the computation of free energy 
differences may be used. We resorted to thermodynamic integration,— in which case 
the function is rewritten as the integral of some canonical averages: 



^(Ai,A2 



— (x, 0)(ix+ / — (Ai,x)dx, 

'^Ai Jo '^A2 



with 



and 



(Ai,A2) = /3 — 

Ai,A2 ) 



dX 



1 + X2ST 



dA 



(Ai,A2) = ((a;-V,.f/((l + Ai)x,2/,z)»^^ 



A2' 



(A2) 



(A3) 
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In conclusion, the numerical procedure consists in first estimating the derivatives of the 
function A and the average potential energy, for as many points as required on the thermo- 
dynamic path chosen. Approximations of S can then be obtained thanks to (lAip . after a 
numerical integration to obtain A. The entropy difference along the path is then plotted, 
and fixing the volume change Ai, the temperature variation is chosen such that the entropy 
difference is 0. This determines A2 as a function of Ai. 



APPENDIX B: RELATIONSHIP BETWEEN THE HUGONIOT AND THE 
ISENTROPE CURVES AT THE POLE 

We present in this Appendix two proofs of the fact that the isentrope curve and the 
Hugoniot agree at order 3 in the volume change. The first one is a standard thermodynamic 
proof, but the second one, based on statistical physics relations, is new to the best of our 
knowledge. 



1. Standard thermodynamic proof 



For the sake of completeness, we reproduce here the proof of Ref. |31|. From the thermody- 
namic relation TdS = dlA+V dV and from the Hugoniot relation^/ = Wo + |('P+'Po)(^o~^); 
the entropy variation along the Hugoniot curve can be computed. One derivation leads to 

A second derivation gives 

' dV^ ' dV dV ^ 2^'" ' ' dV^ 



With a final derivation, 

d^Sung , d'^SuvLgdTuug , dSungd'^Tung _ id'^Vnug , 1... ^.-.d^Vung ,-^^„^ 

^^v^ + ^^v^^y- + ~l\r^v^ - -2~ly^ + 2^^° " ^^^v^- ^^^^ 

At the initial state (denoted with a subscript 0), that is, in the limit V ^ V^/it holds: 



d^S, 



dV 

'Hug 
2 



0, 



dV 



18 



The entropy variation along the Hugoniot curve is therefore of order three in volume, and 
the Hugoniot and the isentropic curves are osculatory. 



2. A statistical physics derivation 

Without loss of generality (and for notational simplicity), we may set 7i(0, 0) = S{0, 0) = 
since we are only interested in differences of JF and S. 

a. Some useful relations. 

The derivatives of the function A are useful for comparing the Hugoniot and the isentrope 
relations. The average xx component of the pressure tensor for the volume V^(Ai) and the 
temperature T(A2) is obtained by averaging the observable 

Therefore, 

-C/(g)/fcBT(A2) 



/ x-V,f/(g)e- 

1 Jv(Xi)'^ 



N 1 + \26T 1 + Ai Jvio) 



'V(Ai)'V 

X ■ V^f/((1 + Ai)a;, y, z) g-^^^i.^aW dq 



(3V{Q) 1 + Ai V{Xi) 









This shows that, using flASp 



(P..)y(A,),T(A.) - + -JVWd^. 

b. Hugoniot curve. 

With the above computations, it is easily seen that the Hugoniot relation (IHl) can be 
restated as 

37V 

/3H(Ai,A2) = —X2ST + (3[{1 + X2ST){{Ux,m))x,m-{{Uo,o))o,o] 

NXifl + X2ST \ ^ifdA.^^. ,^ ^ ^^,dA.^ , A 
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c. Comparison between the Hugoniot and the isentrope. 



We now Taylor expand the difference j3TC(Xi, X2) — 5(Ai, A2) up to the third order, i.e. 
neglecting a remainder term r(Ai,A2) which is such that |r(Ai,A2)| < C'dAil + A2|)'^. We 
denote such remainders by ©(A"^) in the sequel. It holds 



3N 



m{K A2) - 5(Ai, A2) = — (A25T - ln(l + X^ST)) + 1 



Ai 



1 



- Ai A2 — 



1 + Ai 



1 + Ai 

l3X26T{{Ux^,X2))xiM 



ln(l + Ai; 



+y (|^(0,0) + |^(A„A2))-A(Ai,A2). 



Introducing the notation 



^. = ^(0,0), 



A, 



dxA, 



(0,0), 



the Taylor expansions of the function A and its first derivatives at an arbitrary state (Ai, A2) 
read (using A{0,0) = 0): 

A(Ai, A2) = AiA + A2A2 + -^An + A1A2A12 + y^2 + O(A^), 



OA 



(Ai, A2) = Ai + XiAn + X2Ai2 + 0{X'^). 



With these Taylor expansions and the relation ( 1A2I) . it is straightforward to show that 



6T 



OA 



WAi,A2)-5(Ai,A2) = — A^5r^ + AiA2— (iV + Ai)+A2(l + A25r) — (Ai,A2) 

4 z 



XoST 



xAn + Ai 



A 
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6T 



+ A2 



3A^ . ^22 



Using (1A2I) . the derivatives of the entropy differences can be computed: 
dS ^ ^ _ N dA I + X2ST d^A 

dx;^^'' - TTx; + dx;^^'' + ~ij^dx;dX2^^'' 

dS 3N6T OA 1 + A2<5T 

8X2^^'^ - 2(1 + A25r) + + ~l^d^2^^'' 

This shows that 

/3H(Ai, A2) - 5(Ai, A2) = ^ (^1^(0' 0) + ^2||(0, 0)) + 0(A='), (B4) 
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so that, since 

dS dS 
S{Xi, X2) = S{0, 0) + Ai^(0, 0) + A2^(0, 0) + 0(A2) 



and 5(0, 0) = 0, it holds 



/?^(Ai, A2) - ( 1 + ^ ) 5(Ai, A2) = O(A^). (B5) 



2 

This relation shows immediately that 7i(Ai, A2) = O(A^) on the isentrope, and so, the initial 
slopes of the curves, and their first derivatives, coincide. 



APPENDIX C: PRECISIONS ON THE ISENTROPIC INTEGRATION 

Several numerical schemes may be used to integrate ffTTl) . The simplest one consists in 
approximating the integral appearing in the exponential factor with a Riemman formula 
using the value of the integrated function on the left side of the interval: 



{V2 -V,)]. (CI) 



Vi 



Finite differences may be 



Of course, higher order integration methods could be used. 

dV 

It remains to decide how to compute the derivative 7: — 

ou 

used to this end, but this would require at least two very carefully converged simulations 
with volumes Vi ± AV. It seems more appealing to compute the partial derivative using 
standard fluctuations formulas:—!^ 



dU 
df 



and 



Vi 

dV 

df 



Vi,Ti 



Vi,Ti / ' 



(C2) 



(C3) 



where C„(\4;^i) is the specific heat at constant volume, and the pressure observable for a 
simulation domain of volume Vi reads 



1 ^ 2 



3V1 ^ mi 

1=1 



(C4) 



The partial derivative dV/dlA can then be evaluated in a single simulation at (A^, Vi,Ti) 
using (IC2|)-(|C3|). 
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The numerical implementation of this method is done as follows. The partial derivative of 
the pressure with respect to the energy is first computed with a Monte Carlo simulation for 
the given initial conditions (A^, Vi, Ti). The temperature T2 is then evaluated from Eq. (ICip . 
The partial derivative is next computed at volume V2 to predict the next temperature. 
Proceeding incrementally, the whole isentrope curve can be constructed. 

The numerical results presented in this work have been obtained by performing canonical 
samplings with a Metropolis algorithm, using the Monte-Carlo Gibbs code.— Partial deriva- 
tives have been computed in the NVT ensemble. The convergence of simple thermodynamic 
averages was generally obtained after A'^gteps = 10^ iterations, but derivative properties (re- 
lated to the covariance of some observables) required about A'steps = 10^ iterations for a 
satisfactory convergence. Error estimates on the canonical samplings have been obtained 
with block averaging,— and the error propagation estimated along the integration scheme 
has been computed using standard propagation rules. In all the cases considered, the sta- 
tistical error (as measured using the 95% confidence interval associated with the variance 
computed from block-averaging) on the predicted temperature on the isentrope curve is 
inferior to 1.5 %. 
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Since the Hugoniot curve and the isentrope are close enough, we assume that the impact of the 
viscosity has roughly the same amplitude on the isentrope curve. 

The Monte-Carlo Gibbs code is owned by the Institut Francais du Petrole, the Universite Paris- 
Sud, and the CNRS, and developped in collaboration with the CEA. 
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Caption list: 

Figure [D (a) (color online) Non-Equilibrium Molecular Dynamics of isentropic release 
waves. The four pictures represent snapshots of the system during the release process, the 
expansion proceeding in the longitudinal x direction. Atoms are colored according to their 
potential energies (scaling corresponding to —1.38 x 10"^'' J for blue up to 7.55 x 10"^'' J 
for red). 

Figure [U (b) (color online) Density profiles taken at different times of the simulation 
(from blue to red as time increases). 

Figure [2j Path in the (Ai,A2) space used to compute states with the same entropy as 
the initial state. Each cross represents some equilibrium canonical sampling along the 
thermodynamic path. First, the isothermal expansion is performed (horizontal line in the 
diagram), starting from the initial state (0,0), until the required density is reached. The 
entropy of the state (Ai, 0) is S'init + A^cxpansion- Then, an isochore cooling is performed (Ai 
is kept fixed; vertical line in the diagram), until the entropy difference during this process is 
the opposite of the entropy variation found in the expansion part. The final state (Ai, A2), 
located at the intersection of the curve AS* = and the vertical line, has then the same 
entropy as the initial state. 

Figure [3l (color online) Isentropic release in a (P, p) diagram. Symbols represent results 
from equilibrium methods, red diamonds for the successive hugoniostat, blue squares for 
the thermodynamic integration and yellow triangles for the entropy integration. NEMD 
results are plotted in green, the width of the so-obtained tube corresponding to the error 
bars. The arrow indicates the path followed during the release. 

Figure HI (color online) Isentropic release in a (T, P) diagram. The symbols are the same 
as in Figure [31 Notice that there is slight deviation of the NEMD results for the lowest 
temperatures. 

Figure [5l (color online) Isentropic release in a (T, p) diagram. The symbols are the 
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same as in Figure [31 There is a noticeable deviation of the NEMD results for the lowest 
temperatures. 
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